Demonstration of three-dimensional contact point determination and contour reconstruction during active whisking behavior of an awake rat

The rodent vibrissal (whisker) system has been studied for decades as a model of active touch sensing. There are no sensors along the length of a whisker; all sensing occurs at the whisker base. Therefore, a large open question in many neuroscience studies is how an animal could estimate the three-dimensional (3D) location at which a whisker makes contact with an object. In the present work we simulated the shape of a real rat whisker to demonstrate the existence of several unique mappings from triplets of mechanical signals at the whisker base to the three-dimensional whisker-object contact point. We then used high speed video to record whisker deflections as an awake rat whisked against a peg, and used the mechanics resulting from those deflections to extract the contact points along the peg surface. These results demonstrate that measurement of specific mechanical triplets at the base of a biological whisker can enable 3D contact point determination during natural whisking behavior. The approach is viable even though the biological whisker has non-ideal, non-planar curvature, and even given the rat’s real-world choices of whisking parameters. Visual intuition for the quality of the approach is provided in a video that shows the contour of the peg gradually emerging during active whisking behavior.

The 12 tracked whisker shapes all varied slightly from each other due to tracking error, so they were averaged together. To perform the averaging, the fifth-order whisker shapes were extended to the greatest x-distance of the longest whisker, and outlier whiskers were eliminated. The remaining whiskers were averaged in the x-y and x-z views and then merged once again to form the final "average" undeflected whisker shape. In whisker-centered coordinates, the origin is placed at the whisker base, and the x-axis is collinear with the proximal region of the whisker. The proximal 70% of the whisker (which is close to planar) is adjusted to lie as close as possible to the x-y plane. The intrinsic curvature of the whisker points the whisker tip in the positive y-direction. B: The whisker-object ("wobj") contact point location in whiskercentered coordinates is given by the coordinates rwobj, θwobj, and Φwobj. The variable rwobj is the linear distance from the origin to the contact point location. θwobj, is the azimuthal angle of the contact point from the x-axis, and Φwobj is the elevation angle of the contact point measured from the x-y plane.

Background on the whisker bending model, Elastica3D
Forces and moments at the base of the whisker for the full range of contact points were computed using a well-established mechanical model, which we call "Elastica3D" [1][2][3][4], found at https://github.com/SeNSE-lab/DigitalRat. The physical parameters used were based on typical whisker values: 100 μm base diameter and a base-to-tip radius ratio of 15 [6,9,10], yielding a tip diameter of 6.67 μm. Young's modulus was set to 3 GPa [11] and the shear modulus G was obtained by inserting the value of Poisson's ratio for keratin, 0.38 [12], into the following equation: The model, Elastica3D, has already been described in detail in previous studies [1][2][3][4], but we reiterate the main principles here. The model takes as input the undeflected shape of the whisker and a contact point location. It then deflects the whisker by applying a point load of magnitude Fapplied at an arc length Sapplied from the base and at an angle zapplied about the whisker's axis. As described in the following paragraphs, Elastica3D runs an optimization over these three variables until the location of the load on the deflected whisker is coincident with the specified contact point location.
For each iteration of the optimization, Elastica3D applies a point force described by Fapplied, sapplied, and zapplied to the whisker. According to quasistatic limits, the applied force remains perpendicular to the whisker's axis at the point of contact.
This force generates a bending moment, MB, and a torque, MX, at each node, and each node rotates in two directions by amounts dictated by linear elastic beam theory. The node rotates about an axis perpendicular to the distal link's axis by an amount described by the following equation: where ds is the length of the link, E is the Young's modulus, and I is the area moment of inertia calculated by I = πr 4 /4, where r is the radius of the whisker at that node.
The node also rotates about axis defined by the distal link by the amount described in the following equation: where ds is again the length of the link, and J is the polar moment of inertia calculated by J = πr 4 /2, and r is the radius of the whisker at that node.
We also impose a rigid boundary condition at the base, meaning that the basepoint does not translate and the proximal-most link does not rotate.
The optimization algorithm used in Elastica3D is the Nelder-Mead algorithm, provided by the fminsearch function in MATLAB, and the cost function of the optimization is the Euclidean distance from the location of the point force on the deflected whisker to the user-specified contact point location. When this distance becomes zero, Elastica3D has found a solution for the whisker and contact point location and outputs both the deflected shape of the whisker and the resultant forces and moments at the whisker base. If the optimization cannot find a solution where the output of the cost function is zero, then Elastica3D determines that the whisker cannot reach the input contact point and declares that the whisker has "slipped off" the contact.

Methods for determining mapping uniqueness
Mapping uniqueness was determined both through visual inspection and neural network analysis. As depicted in Figure 3 of the main text, the mappings can be visualized as a series of monochromatic surfaces; the colors of the surfaces represent the value of rwobj, θwobj, or Φwobj for a given combination of mechanical variables. If any of these monochromatic surfaces overlap, then the mapping is not unique: the overlap implies that a single reading of the three mechanical variables used would give two or more conflicting results for the contact point location. Thus, in order to pass the visual inspection for uniqueness, a mapping could not have any apparent overlap.
The second uniqueness validation involved using a neural network as a function approximator. In this approach, a neural network was constructed to solve for the highly nonlinear function that related the three selected mechanical variables to the contact point, rwobj, θwobj, and Φwobj. If the neural network could not solve for the nonlinear function within certain limitations (described in detail below), then the mapping was determined to be non-unique.
We used the MATLAB built-in neural network functionality to create, train, and run the neural networks. All networks used the same architecture, selected to solve best for highly nonlinear functions. All networks had three input nodes for the three mechanical signals, three output nodes for the three whisker-object contact point coordinates, and two hidden layers of ten nodes each. The nodes in the hidden layer used a hyperbolic tangent sigmoid function, and the nodes in the output layer used a linear transfer function. The networks were trained using Levenberg-Marquardt backpropagation, and the performance function was the mean squared normalized error. The training set encompassed 35% of the data. MATLAB's built-in training functions used 7.5% of the remaining data for testing during training and a further 7.5% for training validation. The remaining 50% of the data was used to calculate errors used for uniqueness determination.
The neural network error was defined as the difference between the three variables rwobj, θwobj, and Φwobj and those same three variables as output by the neural network. These errors were then weighted by 2 ( ) to determine the median errors in rwobj, θwobj, and Φwobj. This weighting was required to ensure that the contact points, which were equally distributed in spherical coordinates, had equal weight in Cartesian coordinates. Mappings were considered unique only if the median weighted neural network errors were less than or equal to 1.5 mm for rwobj, and less than or equal to 1 degree for θwobj and Φwobj.
We found that some mappings were unique only in select regions of the space. For example, some mappings were unique only when points that involved large deflections were excluded. Contact points were determined to be large deflection if the link that rested on the contact point was oriented in the negative x-direction. Other mappings were unique only when points that resulted from either "concave forward" or "concave backwards" collisions were excluded. Points were determined to be concave forward if, in the top-down view in whisker-centered coordinates, they lay below (in the negative y-direction) the points that defined the curve of the whisker. Points above the whisker (in the positive y-direction) were considered concave forward.
Finally, it is important to note a subtlety in the term "uniqueness." Neural network errorsaveraged across the entire mapping -could be well below threshold for all three geometric coordinates, but still have tiny regions of extreme non-uniqueness. This situation occurred for the case of the mapping used throughout this study (FX, MB, MD), resulting in less accurate rwobj estimates for points very near the whisker, as further discussed in the text for Figure 5 and 6.

Obtaining estimates of mapping error
As described in the previous section, estimates of mapping uniqueness were based on computing the "neural network error." The neural network error is a measure of the performance of the neural network, i.e., its ability to predict 50% of the data (the test set) based on the other 50%.
However, errors between actual contact point and predicted contact point are present throughout the entire mapping, and thus a different metric is required to quantify the error and resolution of the mapping. The "mapping error" is found using 100% of non-normalized data.
To find the mapping error, all combinations of FX, MB, and MD from the entire range of whiskerobject contact points were input into the neural network, which calculated contact point locations based on these mechanical signals. The contact point locations estimated by the neural network were then compared to the actual contact point locations, and the Euclidean distance between these two defined the mapping error at each contact point. In future work, this error could potentially be reduced by optimizing the neural network architecture and associated training algorithm, or perhaps by using a larger training dataset. An alternative possibility is that the large error in the region of small FX and MB might actually reflect nonuniqueness at the resolution analyzed. Although previous work has demonstrated that a perfectly parabolic curve produces a unique mapping with the triplet Fx, MB, and MD, the real whisker used in the present work curves out of the x-y plane, and has irregular curvature.

Mappings using other mechanical variable triplets
We also evaluated uniqueness for the remaining 19 mappings resulting from all the different triplet combinations of the six mechanical variables at the whisker base. The uniqueness results are listed in Table A (next page). Each mapping triplet is listed in the columns on the left, and the rightmost column describes the region of the contact point space in which the mapping is unique. Mappings that are labeled "All" are unique when all contact points are included. Note that the Fx, MB, MD mapping featured in the present work is one of three possible combinations that yields a unique mapping for all contact points.
Some mappings are labeled "ELD," indicating that the mappings are unique only when large deflections are excluded. Large deflections are defined as those that deflect the whisker past the point where FD experiences a 180° flip.
Mappings can also be labeled "CF" or "CB," meaning that the mappings are unique only when either concave forward or concave backward collisions are considered, respectively.
A mapping labeled either "ELD CF" or "ELD CB" means that the mapping is unique when only concave forward or concave backwards collisions are considered and large deflections are excluded.
Only two mappings, on rows 19 and 20, were not unique under any of these conditions. Some mappings are unique under multiple conditions, such as the FT, FD, MD mapping on row 6. It is unique either if large deflections are excluded or if only concave backward collisions are included. To represent this, the right-hand column in Table A (next page) is split into multiple boxes to show both "ELD" and "CB." Mappings can also be marked with a single asterisk. This asterisk indicates that although the mapping is unique for the regions indicated, it is non-unique in the x-y plane. The mapping cannot discriminate between points in this region because MX is always equal or nearly equal to 0.
Although the bulk of this study has focused on the FX, MB, MD mapping, this table indicates alternative mapping combinations that can be further explored in future studies.  Figure 3. A tabulation of the uniqueness for all 20 mappings. Boxes marked "All" indicate that the mapping is unique for the entire contact point space. Boxes marked "ELD" indicate that the mapping is unique only after excluding large deflections. Boxes marked "CF or "CB" indicate that the mapping is unique when only concave forward (CF) or concave backwards (CB) contact points are included. A mapping labeled "ELD CF" or "ELD CB" means that the mapping is unique only when concave forward or concave backwards are considered and large deflections are excluded. Some boxes include a forward slash and give multiple readings. This notation indicates that the mapping was unique for multiple conditions. For example, the way to read mapping 12 is that the mapping is unique if it is known whether the points are concave backwards or concave forward, provided large deflections are excluded. The way to read mapping 6 is that it is unique either if large deflections are excluded, or if only concave backwards points are considered. The mappings with boxes marked "not unique" had no well-defined regions of uniqueness. Two of the mappings are marked with an asterisk, indicating that they are non-unique specifically when the contact points are in the x-y plane, where MX is always 0.

Inputs
When is the mapping unique? 1 2 3